Rossmann Sales Forcasting¶
Introduction¶
Rossmann operates over 3,000 drug stores in 7 European countries. Currently, Rossmann store managers are tasked with predicting their daily sales for up to six weeks in advance. Store sales are influenced by many factors, including promotions, competition, school and state holidays, seasonality, and locality. With thousands of individual managers predicting sales based on their unique circumstances, the accuracy of results can be quite varied. You are provided with historical sales data for 1,115 Rossmann stores.
Data Description¶
- Id - an Id that represents a (Store, Date) duple within the test set
- Store - a unique Id for each store
- Sales - the turnover for any given day (this is what you are predicting)
- Customers - the number of customers on a given day
- Open - an indicator for whether the store was open: 0 = closed, 1 = open
- StateHoliday - indicates a state holiday. Normally all stores, with few exceptions, are closed on state holidays. Note that all schools are closed on public holidays and weekends. a = public holiday, b = Easter holiday, c = Christmas, 0 = None
- SchoolHoliday - indicates if the (Store, Date) was affected by the closure of public schools
- StoreType - differentiates between 4 different store models: a, b, c, d
- Assortment - describes an assortment level: a = basic, b = extra, c = extended
- CompetitionDistance - distance in meters to the nearest competitor store
- CompetitionOpenSince[Month/Year] - gives the approximate year and month of the time the nearest competitor was opened
- Promo - indicates whether a store is running a promo on that day
- Promo2 - Promo2 is a continuing and consecutive promotion for some stores: 0 = store is not participating, 1 = store is participating
- Promo2Since[Year/Week] - describes the year and calendar week when the store started participating in Promo2
- PromoInterval - describes the consecutive intervals Promo2 is started, naming the months the promotion is started anew. E.g. "Feb,May,Aug,Nov" means each round starts in February, May, August, November of any given year for that store
Exploratory Data Analysis¶
In [1]:
# Importing libraries
import numpy as np
import pandas as pd
from pandas import datetime
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
import seaborn as sns
import matplotlib.gridspec as gridspec
import plotly.express as px
from IPython.display import display
%matplotlib inline
from sklearn.metrics import mean_squared_error
from statsmodels.distributions.empirical_distribution import ECDF
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.arima_model import ARMA,ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
from pmdarima import auto_arima
from fbprophet import Prophet
import warnings
warnings.filterwarnings("ignore")
In [2]:
# Importing data
train = pd.read_csv("train.csv")
store = pd.read_csv("store.csv")
In [3]:
train.head()
Out[3]:
In [4]:
train.info()
In [5]:
# Change Date datatype to datetime type
train['Date'] = pd.to_datetime(train['Date'])
In [6]:
# Date extraction
train['Year'] = train['Date'].dt.year
train['Month'] = train['Date'].dt.month
train['Day'] = train['Date'].dt.day
train['WeekOfYear'] = train['Date'].dt.weekofyear
In [7]:
# Adding new variable
train['SalePerCustomer'] = train['Sales']/train['Customers']
train['SalePerCustomer'].describe()
Out[7]:
On average customers spend about 9.50$ per day.
In [8]:
# ECDF: empirical cumulative distribution
sns.set(style = "ticks")
plt.figure(figsize = (12, 6))
plt.subplot(311)
cdf = ECDF(train['Sales'])
plt.plot(cdf.x, cdf.y, label = "statmodels");
plt.title('Sales'); plt.ylabel('ECDF');
# plot second ECDF
plt.subplot(312)
cdf = ECDF(train['Customers'])
plt.plot(cdf.x, cdf.y, label = "statmodels");
plt.title('Customers');
# plot second ECDF
plt.subplot(313)
cdf = ECDF(train['SalePerCustomer'])
plt.plot(cdf.x, cdf.y, label = "statmodels");
plt.title('Sale per Customer');
plt.subplots_adjust(hspace = 0.8)
About 20% of data has zero amount of sales/customers probably due to the fact that the store is closed for state holidays or school holidays
Data Wrangling¶
In [9]:
# Copy data for cleaning
df_sales = train.copy()
df_store = store.copy()
Train data¶
In [10]:
# Closed Stores with zero sales
df_sales[(df_sales.Open == 0) & (df_sales.Sales == 0)].head()
Out[10]:
In [11]:
len(df_sales[(df_sales.Open == 0) & (df_sales.Sales == 0)])
Out[11]:
There're 172817 closed stores in the data. It is about 10% of the total amount of observations.
In [12]:
# Opened stores with zero sales
df_sales[(df_sales.Open == 1) & (df_sales.Sales == 0)].head()
Out[12]:
In [13]:
len(df_sales[(df_sales.Open == 1) & (df_sales.Sales == 0)])
Out[13]:
There're 54 opend stores but no sales in the data.To avoid any biased forecasts we will drop these values.
In [14]:
# Drop data
df_sales = df_sales[(df_sales["Open"] != 0) & (df_sales['Sales'] != 0)]
In [15]:
df_sales.head()
Out[15]:
In [16]:
df_sales.shape
Out[16]:
Store data¶
In [17]:
df_store.head()
Out[17]:
In [18]:
df_store.info()
In [19]:
# Null value
df_store.isnull().sum()
Out[19]:
In [20]:
# Missing values in CompetitionDistance
df_store[pd.isnull(df_store.CompetitionDistance)]
Out[20]:
In [21]:
# Fill NaN with a median value (skewed distribuion)
df_store['CompetitionDistance'].fillna(df_store['CompetitionDistance'].median(), inplace = True)
In [22]:
# Fill other NA by 0
df_store.fillna(0, inplace = True)
In [23]:
# Join sales and store data
df_clean = pd.merge(df_sales, df_store, how = 'inner', on = 'Store')
df_clean.head()
Out[23]:
Data analysis¶
In [24]:
# Sales distribution among stores
df_clean.groupby('StoreType')['Sales'].describe()
Out[24]:
StoreType B has the highest average of Sales among all others, however we have much less data for it.
In [25]:
# Total sales and customers among stores
df_clean.groupby('StoreType')['Customers', 'Sales'].sum()
Out[25]:
StoreType A comes as first, StoreType D goes on the second place in both Sale and Customers
In [26]:
# sales trends
sns.factorplot(data = df_clean, x = 'Month', y = "Sales",
col = 'StoreType',
palette = 'plasma',
hue = 'StoreType',
row = 'Promo')
Out[26]:
In [27]:
# sales trends
sns.factorplot(data = df_clean, x = 'Month', y = "Customers",
col = 'StoreType',
palette = 'plasma',
hue = 'StoreType',
row = 'Promo')
Out[27]:
In [28]:
# sale per customer trends
sns.factorplot(data = df_clean, x = 'Month', y = "SalePerCustomer",
col = 'StoreType',
palette = 'plasma',
hue = 'StoreType',
row = 'Promo')
Out[28]:
The plots above showed StoreType B as the most selling store but it is not true. The highest SalePerCustomer amount is observed at the StoreType D, about 12€ with Promo and 10€ without. As for StoreType A and C it is about 9€.
Conclusion of EDA¶
- The most selling and crowded StoreType is A.
- Best "Sale per Customer" StoreType D.
- Low SalePerCustomer amount for StoreType B.
- Customers tends to buy more on Modays when there's one promotion (Promo) and on Sundays when there's no promotion at all (both Promo and Promo1 are equal to 0).
- Promo2 alone doesn't seem to have any significant effect in the Sales.
Time Series Analysis¶
- Time series seems an alternative approach when dealing with continuous values, apart from ML regresions
- Time in time series is dependant while target values is independant.
- Seasonality, Trends, Residuals are 3 major things worth being explored in seasonal_decompose. Plus, some special event, holidays, ect, are captured as well.
- Plotting ACF and PACF is the essential step before choosing which time series model is applyied, such as AR, MA, ARMA, ARIMA, SARIMA.
- Prophet Facebook function is quite easy to use and is similar to ARIMA family.
Seasonality¶
We chose only one strore, Store 2 for StoreType A to demonstrate the analysis It also makes sense to downsample the data from days to weeks using the resample method to see the present trends more clearly.
In [29]:
# Choose store 2
store_2 = df_sales[df_sales['Store'] == 2][['Sales', 'Date']]
In [30]:
fig = px.line(store_2,x='Date',y ='Sales')
fig.update_xaxes(
rangeslider_visible=True,
rangeselector = dict(
buttons=list([
dict(count=1, label='1y',step='year',stepmode='backward'),
dict(count=2, label='2y',step='year',stepmode='backward'),
dict(count=3, label='3y',step='year',stepmode='backward'),
dict(step='all')
])))
fig.show()
In [31]:
store_2['Date'].sort_index(ascending = False, inplace=True)
In [32]:
store_2 = store_2.set_index('Date').resample('W').mean()
store_2.head()
Out[32]:
In [33]:
plt.figure(figsize=(12, 5))
sns.lineplot(x=store_2.index, y=store_2['Sales'])
Out[33]:
Trend, Seasonal, Residual¶
In [34]:
f , (ax1, ax2, ax3, ax4)= plt.subplots(4, figsize = (12, 10))
decomposition = seasonal_decompose(store_2['Sales'], model = 'multiplicative')
decomposition.observed.plot(ax = ax1)
ax1.set_ylabel('OBSERVED')
decomposition.trend.plot(ax = ax2)
ax2.set_ylabel('TREND')
decomposition.seasonal.plot(ax = ax3)
ax3.set_ylabel('SEASONAL')
decomposition.resid.plot(ax = ax4)
ax4.set_ylabel('RESIDUAL')
f.subplots_adjust(hspace = 0.5)
Downward trend with seasonal effect
Autocorrelation and Partial Autocorrelation Plot¶
In [35]:
fig = plt.figure(constrained_layout=True, figsize=(15, 4))
grid = gridspec.GridSpec(nrows=1, ncols=2, figure=fig)
ax1 = fig.add_subplot(grid[0, 0])
plot_acf(store_2['Sales'], lags = 20, ax=ax1);
ax2 = fig.add_subplot(grid[0, 1])
plot_pacf(store_2['Sales'], lags = 20, ax=ax2);
plt.show();
- ACF is a measure of the correlation between the timeseries with a lagged version of itself. I choose lags of 20. ACF help us choose MA model's parameter.
- PACF, on the other hand, measures the correlation between the timeseries with a lagged version of itself, after removing the influence of any variance in the middle. PACF helps us choose AR model's parameter.
- Looking at the Autocorrelation plot shows store 2 has stationarity in sales data.
Time Series Modeling¶
In [36]:
# Double check data stationary by using adfuller method.
from statsmodels.tsa.stattools import adfuller
def adf_test(series,title=''):
"""
Pass in a time series and an optional title, returns an ADF report
"""
print(f'Augmented Dickey-Fuller Test: {title}')
result = adfuller(series.dropna(),autolag='AIC')
labels = ['ADF test statistic','p-value','# lags used','# observations']
out = pd.Series(result[0:4],index=labels)
for key,val in result[4].items():
out[f'critical value ({key})']=val
print(out.to_string())
if result[1] <= 0.05:
print("Strong evidence against the null hypothesis")
print("Reject the null hypothesis")
print("Data is stationary")
else:
print("Weak evidence against the null hypothesis")
print("Fail to reject the null hypothesis")
print("Data is non-stationary")
In [37]:
adf_test(store_2['Sales'])
The data is good enough to perform time series models.
In [38]:
# Process test data
train = store_2[:120]
test = store_2[120:]
train.head()
Out[38]:
ARMA model¶
In [39]:
# auto_arima help us choose the optimal model, sometime manually tweaking model hyperparameters yeild better result.
model=auto_arima(train['Sales'],m=7,seasonal=False,
start_p=0, start_q=0, max_order = 5, test = 'adf',
error_action='ignore', suppress_warning=True,
stepwise=True, trace = True)
In [40]:
model.summary()
Out[40]:
In [41]:
#ARMA
model1 = ARMA(train['Sales'], order=(2,0))
results=model1.fit()
start=len(train)
end=len(train)+len(test)-1
predictions1 = results.predict(start=start, end=end)
predictions1.index = test.index
test['Sales'].plot()
predictions1.plot(label='prediction');
print('Root Mean Squared Error: ', np.sqrt(mean_squared_error(test['Sales'], predictions1)))
SARIMA model¶
In [42]:
model=auto_arima(train['Sales'],m=7,seasonal=True,
start_p=0, start_q=0, max_order = 5, test = 'adf',
error_action='ignore', suppress_warning=True,
stepwise=True, trace = True)
In [43]:
model.summary()
Out[43]:
In [44]:
model = SARIMAX(train['Sales'],order=(2, 0, 0),seasonal_order=(0, 0, 1, 7))
results = model.fit()
start=len(train)
end=len(train)+len(test)-1
predictions1 = results.predict(start=start, end=end, dynamic=False)
predictions1.index = test.index
test['Sales'].plot()
predictions1.plot(label='prediction');
plt.legend()
print('Root Mean Squared Error: ', np.sqrt(mean_squared_error(test['Sales'], predictions1)))
FB Prophet¶
In [45]:
df_prophet = store_2.reset_index()[['Date','Sales']].rename({'Date':'ds','Sales':'y'}, axis='columns')
df_prophet.head()
Out[45]:
In [46]:
df_prophet.tail(20)
Out[46]:
In [47]:
train_prophet = df_prophet[:120]
test_prophet = df_prophet[120:]
In [48]:
test_prophet.head()
Out[48]:
In [58]:
model = Prophet(interval_width=0.95,yearly_seasonality=True,weekly_seasonality=True)
In [59]:
model.fit(train_prophet)
Out[59]:
In [60]:
future = model.make_future_dataframe(periods=15,freq='W')
future.tail()
Out[60]:
In [61]:
forecast = model.predict(future)
forecast[['ds','yhat','yhat_lower','yhat_upper']].tail(20)
Out[61]:
In [62]:
pd.concat([df_prophet.set_index('ds')['y'],forecast.set_index('ds')['yhat']],axis=1).plot()
Out[62]:
In [63]:
model.plot(forecast);
In [64]:
model.plot_components(forecast);
In [57]:
from fbprophet.plot import add_changepoints_to_plot
fig = model.plot(forecast)
a = add_changepoints_to_plot(fig.gca(),model,forecast)
LSTM¶
In [70]:
df=df_sales[df_sales['Store'] == 2][['Sales', 'Date']]
df['Date'].sort_index(ascending = False, inplace=True)
In [71]:
from sklearn.preprocessing import MinMaxScaler
scaler=MinMaxScaler(feature_range=(0,1))
data=scaler.fit_transform(df[['Sales']])
In [80]:
train = int(len(data)*0.75)
train_data,test_data=data[:train],data[train:]
In [85]:
# convert an array of values into a dataset matrix
def create_dataset(dataset, look_back=1):
dataX, dataY = [], []
for i in range(len(dataset)-look_back-1):
a = dataset[i:(i+look_back), 0]
dataX.append(a)
dataY.append(dataset[i + look_back, 0])
return np.array(dataX), np.array(dataY)
In [86]:
# reshape into X=t and y=t+1
look_back = 1
train_X, train_y = create_dataset(train_data, look_back)
test_X, test_y = create_dataset(test_data, look_back)
In [87]:
train_X.shape, train_y.shape
Out[87]:
In [88]:
# reshape input to be [sample, time steps, features] which is required for LSTM
train_X =train_X.reshape(train_X.shape[0], 1, train_X.shape[1])
test_X = test_X.reshape(test_X.shape[0], 1, test_X.shape[1])
In [93]:
### Create the LSTM model
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense
from tensorflow.keras.layers import LSTM
model=Sequential()
model.add(LSTM(50,return_sequences=True,input_shape=(1,1)))
model.add(LSTM(50))
model.add(Dense(1))
model.compile(loss='mean_squared_error',optimizer='adam')
model.fit(train_X,train_y,validation_data=(test_X,test_y),epochs=100,batch_size=1,verbose=1)
Out[93]:
In [94]:
#Model Prediction
train_predict = model.predict(train_X)
test_predict = model.predict(test_X)
In [95]:
#Transforming data back to original form
train_predict=scaler.inverse_transform(train_predict)
test_predict=scaler.inverse_transform(test_predict)
In [97]:
## Calculate RMSE performance metrics
from sklearn.metrics import mean_squared_error
import math
math.sqrt(mean_squared_error(train_y,train_predict))
Out[97]:
In [ ]: